Hydrodynamic Shock Wave Studies within a Kinetic Monte Carlo Approach 



Irina Sagert 1 , Dirk Colbry 2 , Terrance Strother 3 , Rodney Pickett 2 , Wolfgang Bauer 1,2 
1 Department of Physics and Astronomy, Michigan State University, East Lansing, Michigan, 4-8824, USA 
^Institute for Cyber- Enabled Research, Michigan State University East Lansing, Michigan 48824, USA 
3 XTD-6, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA 

(Dated: October 31, 2012) 

Kinetic approaches are routinely employed to simulate the dynamics of systems that are too 
rarified to be described by the Navier-Stokes equations. However, generally they are far too compu- 
tationally expensive to be applied for systems that are governed by continuum hydrodynamics. In 
this paper, we introduce a massively parallelized test-particle based kinetic Monte Carlo code that 
is capable of modeling the phase space evolution of an arbitrarily sized system that is free to move 
in and out of the continuum limit. Using particle mean free paths which are small with respect to 
the characteristic length scale of the simulated system, we retrieve continuum behavior, while non- 
equilibrium effects are observed when the mean free path is increased. To demonstrate the ability of 
our code to reproduce hydrodynamic solutions, we apply a test-suite of classic hydrodynamic shock 
problems. Simulations using tens of millions of test-particles are found to reproduce the analytical 
solutions well. 

PACS numbers: 47.45.Ab,45.50.-j,47.85.Dh, *43.25.Cb,07.05.Tp, 82.20.Wt 

I. INTRODUCTION AND MOTIVATION 

In fluid dynamics, when the characteristic length scale of a flow does not greatly exceed the mean free path of the 
particles it is comprised of, the continuum approximation assumed by the Navier-Stokes (NS) equations breaks down 
and the particle nature of matter must be explicitly taken into account pQ. Under these conditions, flows are said to 
be rarified. The rarefaction of a flow is characterized by the Knudsen number Kn, which is defined as the ratio of the 
mean free path A to a characteristic length scale of the system L: 

Kn = X/L. (1) 

While the continuum limit of hydrodynamics is applicable when Kn -C 1, it is generally accepted that flows with 
Kn > 0.01 are not well described by the NS equations as they do not form a closed set in this regime [Tj . There are cur- 
rently many areas of research in which flows with large Knudsen numbers are relevant. The studies of hypersonic flow 
[5] and shock structure [TJ, simulation studies of space flight [3], nano-scale devices [I], particle production in heavy 
ion collisions [5]-[7], the dynamics of inertial confinement fusion (ICF) capsules [BlU0|. as well as out-of-equilibrium 
neutrino-matter interactions in core-collapse supernovae [11] all require the ability to model flows with Kn > 0.01. 
The latter two subjects are of particular interest to us. 

Despite the enormous efforts poured into achieving thermonuclear ignition of ICF capsules at the National Ignition 
Facility (NIF), satisfactory yield has never been obtained. The ICF capsule experiments at NIF are guided by complex 
numerical simulations. All of the codes that suggest that the techniques employed thus far to ignite the nuclear fuel 
should have been successful use continuum hydrodynamics to describe the entire capsule system. The so-called "hot 
spot" of a typical ICF capsule, where thermonuclear burn is expected to initiate and shock fronts must be resolved, is 
approximately 25 micrometers across and has a density and temperature on the order of 100 g/cc and 10 keV respec- 
tively [12 . It is a trivial exercise to demonstrate that in this environment, the mean free path of a thermal deuteron, 
a critical component of the nuclear fuel, is on the order of a few micrometers. Clearly the continuum approximation 
is not applicable here. Furthermore, investigations of shock front kinetic effects depleting the ion distribution of fast 
particles in ICF capsules [5 HTU1 fT3HTrj] found that the resultant highly non-Maxwellian ion distributions may alter 
the hydrodynamics and thermonuclear yield of the capsules. Other studies suggest that capsule yield can also be 
diminished by kinetic effects associated with strong self-generated electric fields which can broaden the shock fronts 
[PTl [T5] . Consequently, there is a strong motivation to develop a kinetic model capable of fully resolving shock fronts 
in all regions of ICF capsules. 

Many theories have been suggested to explain the core-collapse supernova explosion mechanism [T9H24] . The determi- 
nation of the role that neutrino-matter interactions play in the generation of these explosions remains an outstanding 
problem in physics (see e.g. [S3] and references therein). Currently, the most commonly accepted model for core- 
collapse supernova explosions is based on deposition of energy by neutrinos that are produced in the central high 
temperature and density region of the collapsed stellar core into cooler less dense regions of the infalling matter 
[T§1 |2"U] . The simulations suggest that this so-called neutrino heating can power fluid instabilities that grow to large 
amplitudes and eventually power an explosion [261 127] . To study such a scenario, an accurate description of the 
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supernova hydrodynamics as well as modeling neutrino transport in three dimensions seems to be ultimately required 
[28H31j . For one-dimensional simulations, it is possible to solve the full Boltzmann equations for neutrino propagation 
J52H31]- However, in two or three dimensions, this approach becomes computationally too expensive and approxima- 
tions have to be made [35 39 . Monte Carlo neutrino transport has been suggested as an alternative approach since it 
is expected to scale better for multi-dimensional simulations (see [40] and references therein). This is one motivation 
for us to develop a large scale Monte Carlo transport code which can be used in supernova studies. Furthermore, most 
simulations operate with nuclear equations of state which evolve only one representative heavy nucleus instead of a full 
statistical ensemble of present nuclei. The inclusion of the latter could alter weak interaction processes and thereby 
the efficiency of neutrino heating as well as impact supernova nucleosynthesis [23 E- 43 . It has been demonstrated 
in previous works that a kinetic approach has the capability to treat baryon and neutrino dynamics identically and 
thereby evolve a full ensemble of nuclei and readily model out-of-equilibrium neutrino-matter interactions |1U I44j . 
Preliminary studies conducted with such an approach found that out-of-equilibrium neutrino-matter interactions with 
nuclei close to the neutron drip line could significantly impact the dynamics of the core-collapse and subsequent ex- 
plosion [TT]|44]. More detailed calculations are necessary to confirm the observed effects. 

It is the long-term goal of this work to contribute to the modeling of systems which are influenced by out-of-equilibrium 
phenomena, such as ICF capsules and core-collapse supernovae. However, for realistic simulations, it is necessary to 
develop an approach capable of resolving macroscopic three-dimensional hydrodynamic shock fronts and describ- 
ing non-continuum flows. The latter requirement necessitates abandoning the NS equations and adopting a kinetic 
description. 



II. THE BOLTZMANN EQUATION 



For non-continuum flows, the governing equation is the Boltzmann equation of kinetic theory. It is a nonlinear 
integro-differential equation that describes the time evolution of the statistical distributions of particles in a fluid that 
undergo binary collisions. The Boltzmann equation is valid for all Knudsen numbers and, in the limit of small mean 
free paths, leads to the the continuum description [35]. For a fluid comprised of N different species of particles, the 
equation satisfied by the probability distribution function (PDF) of the ith species is [33] : 

9fi(r,p,t) - dfii^ftt) 

1 Vfi(r,p,t)+F — =I»,cou 5 (2) 

at m op 

where r, p, and mi are the position, momentum, and mass of particle i respectively, F is the external force field, 
fi(f,p, t) is the number of particles of species i found in a differential phase space neighborhood of f and p at time t, 
and l^coii is the term that takes into account the changes in fa induced by two-body collisions. ij )CO ii is generally a 
sum of integrals involving the pre-collision and post-collision PDFs of all species present |3S] . It is this collision term 
that couples the N equations and generally makes them very difficult to solve. 



A. Analytic Solution techniques 



While the Boltzmann equation in not generally solvable by analytical techniques, approximations can be employed 
to derive less complex equations. By making assumptions about the form of the PDFs, the Boltzmann equation can 
be used to derive approximate fluid dynamics equations of higher order in Knudsen number than the NS equations [5] . 
These higher-order fluid dynamic models are the so-called extended or generalized hydrodynamics equations. Another 
common approximate technique is to modify the Boltzmann equation itself. The most famous example of this is the 
Bhatnagar-Gross-Krook (BGK) [57] equation that replaces the collision term in eq.([2]) with a relaxation term: 

T f~M,i ft , \ 
^i.coll — r \o) 

r 

where fu,i is the Maxwell distribution for the ith species, fi is the PDF of the ith species, and r is a particle collision 
time. Despite the fact that the BKG equations are less complicated than the full Boltzmann equation, numerical 
techniques such as the lattice Boltzmann method or discrete velocity method [47H50] are typically required to solve 
them. 
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B. Numerical Simulation Techniques 



In the absence of a general analytic solution technique for the full Boltzmann equation, simulation schemes that 
attempt to directly model the dynamics of a fluid using a large number of simulated particles have been developed. 
While it is possible to infer the PDFs from the simulated particles phase space coordinates, the problem can be 
formulated without reference to them. The first of these direct simulation approaches was the molecular dynamics 
(MD) method [51]. In these calculations, each physical particle is represented by a simulated particle. Many advances 
have been made in MD calculations providing unique insight to fluid dynamics on molecular scales [52-54 . However, 
the one-to-one physical-to-simulated particle ratio computationally limits MD calculations to volumes much smaller 
than many hydrodynamic features of interest. 



C. The Test-Particle Method 



In order to model the phase space evolution of systems that contain too many physical particles to make use of a 
MD-like direct simulation scheme, so-called test-particle methods are employed. These allow the physical-to-simulated 
particle ratio to greatly exceed one and thereby permit the modeling of macroscopic samples of fluids. Rather than 
attempting to discretize the entire phase space relevant to the simulated system and track the value of the phase 
space densities in each cell in time, the test-particle method only tracks the initially occupied phase space cells and 
represents them by imaginary test-particles. These are propagated in a way that models the physical evolution of the 
phase space. The test-particles interact with one another via one-body mean field forces and scatter with realistic 
cross sections. 

Formally, the test-particle method approximates the phase space density with a sum over delta functions [55] 

AT tp 

f(r,p,t) = 5> 3 (r-- m)6 3 (p-Pi(t)), (4) 

i=0 

where -/V tp is the total number of test-particles. The initial coordinates of these delta function test-particles are 
determined by the initial conditions of the simulated system. Inserting this solution into the system transport 
equations such as the Boltzmann equations in eq.([2| generates the following set of simple first-order linear differential 
equations that govern the motion of the i th test-particle's centroid coordinates in the full six-dimensional phase space: 

d -& = F{r t )+C{p t ) 



dt 
d_ 
dt 



d _ fi 

Ti = 



m i 



(5) 



Here F(fi) is the mean field one-body force acting on the i th test-particle, and C(pi) symbolizes the effects that 
two-body collisions with other test-particles have on the i th test-particle's momentum. 

For a system comprised of A^ p h ys physical particles being modeled with N tp test-particles, the ratio -/Vphys/iVtp ef- 
fectively determines a cutoff scale below which details cannot be resolved. Naturally, the test-particle approach is 
applicable as long as N tp is sufficiently large to capture the gross dynamics of the system's phase space. When 
-^Vphys/^Vtp becomes too large, some "small-scale" phenomena can become impossible to resolve. It must therefore 
be established that these details do not impact the gross phase space dynamics and/or can be taken into account 
indirectly. The former can be accomplished with convergence tests. 

However, multiscale problems are certainly not unique to the test-particle approach. Hydrodynamic calculations have 
to spatially discretize simulated systems. When the corresponding volumes are too large the same scale and resolu- 
tion concerns can occur. Representative particle models, which are very similar to test-particle calculations in the 
large -/V p hys/-^tp limit, have to contend with scale and resolution issues as well when the importance sampling of the 
particles the system is comprised of is made. Failure to represent physical particles with certain characteristics with 
a sufficient number of representative particles can prevent calculations from resolving details that are essential to the 
gross system behavior. 



III. THE DSMC TECHNIQUE 



The test-particle method employs probabilistic procedures to model two-body collisions. It is therefore classified as 
a Direct Simulation Monte Carlo (DSMC) technique [TH [SU [57| . primary approximation of the DSMC approach 
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is the decoupling of the motion and the collisions of simulated particles during sufficiently small time intervals. DSMC 
techniques employ an operator splitting method to separate the Boltzmann equation into two processes: transport 
and collisions. Simulated particle motions are calculated deterministically and two-body collisions between them are 
modeled probabilistically. Though variations exist between different implementations of DSMC algorithms, a given 
iteration can be described by four basic processes: sample flow held, move the particles accordingly, organize the 
particles into a scattering grid, and model two-body collisions. The first two processes involve evaluating the flow 
field forces acting on the simulated particles and then solving the resultant first-order linear differential equations of 
motion. The third process consists of dividing the simulated volume into grid cells and organizing the particles into 
them. The size chosen for the scattering grid cells can depend strongly upon the way the two-body collisions are 
modeled in the forth step of the iteration. 

A common approach to implement two-body particle scattering is to model collisions only between particles in the same 
scattering grid cell wherein the collision partners are selected randomly []3J [S5J [SS] • It has been demonstrated that, in 
the limits of infinite simulated particle count and vanishing scattering grid cell and time step sizes, DSMC algorithms 
that employ this type of scattering technique yield results that converge to known solutions of the Boltzmann equation 
[60j . In particular, this scattering technique has been used extensively and with great success in the simulation of 
microscopic flows [551 . 

To avoid unphysical collisions between simulated particles that are too far apart from each other in a given time 
step, the size of the scattering grid cells has to be tightly connected to mean free path of the physical particles. This 
requirement can make scattering algorithms of this type prohibitively expensive for flows with mean free paths and 
particle collision times much smaller than the characteristic scales for lengths and times. For such flows, a way of 
modeling two-body collisions that selects scattering pairs more judiciously is required. 

In this work we apply such an altered DSMC algorithm for the detection and performance of particle collisions. Our 
approach will be described in sections |IV| and |Vj while in section |VI| we present our results for shock wave modeling 



with the kinetic approach. The paper is closed with the conclusions and outlook in sections VIII and IX respectively. 



IV. THE DISTANCE-OF-CLOSEST-APPROACH METHOD FOR SCATTERING PARTNER SEARCH 

Unlike in traditional DSMC algorithms where scattering partners are chosen randomly from a grid cell, in our 
approach the possible collision between two test-particles is determined from their distance of closest approach d m i n . 
If the latter is reached during a timestep At the paths of two particles A and B can cross and a collision is possible. 
To detect an intersection of particle paths, the relative position vectors r re i at the current as well as the next timestep: 

r re i (*) = Pa (t) ~ r B (t) , r rel (t + At) = f A (t + At ) - r B (t + At ) (6) 

have to be projected on the corresponding relative velocity vectors v le \. This results in the crossing number: 

X = (Pvdt) ■ v iel (t))(P Tel {t + At) ■ v Icl (t + At)). (7) 

A negative value of \ indicates that within the time step Ai, the particles have reached d m { n and their paths have 
crossed. This is a first indication of a possible collision. However, its occurrence depends on the interaction cross- 
section cr, i.e. the particle mean free path A. Both can be related to an effective particle interaction radius r e ff 
via: 

A = (cm)" 1 = (4tt r 2 cS n)" 1 , (8) 

whereas n = N/V is the number of particles N per volume V. Two possible scattering partners A and B can 
only interact with each other if their distance of closest approach is within the sum of the effective radii d m j n < 
{f e s,A + ^eff.s)- The collision time t c can then be calculated from: 

|r re i + Wrci t c \ = r e $ tA + r e K,B, (9) 
|wroi| 2 t\ + 2 (v Ici ■ rUi) t c + \f Icl \ 2 = (r cffij4 + r cffjB ) 2 (10) 

to: 




If the collision times t c \ and t c ^ are real numbers a collision can take place. In case the particle distance is larger than 
the sum of the effective radii, that is |r r ci| 3* (feff,^ +r e ff,_B), the particle cross-sections are too small and the collision 
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[100] [200] [300] [400] Parallel neighborhood search by six CPUs 



FIG. 1: (a) Visualization of a 2D particle setup with 25 bins. The active bin is located at the center (light grey), surrounded 
by its 8 neighborhood bins (dark grey), (b) Parallel neighborhood binary collision detection performed by six processors in 
the active bins. The neighboring bins are numbered according to their active bins, (c) Two possibilities for the progress of 
a parallel scattering partner search for six CPUs. (1) Simultaneous update of active bins for all CPUs. (2) Individual CPU 
update of the active bin. The latter can not be part of the neighborhood in the ongoing collision detection of the remaining 5 
processors. 



times are imaginary. Negative values of t c i a, on the other hand, indicate that the effective radii of particles A and B 
are overlapping already at the current time step. This situation can occur frequently for high particle densities or if 
the value of A is chosen much smaller than the distance a test-particle can travel within At. As will be described in 
the next section, we distinguish between positive and negative collision times by choosing a slightly different approach 
for position update of the corresponding particles. 

Typically, the analysis of t c finds several potential collision partners for each particle. Once the latter have been 
determined, the final partners can be chosen either according to the shortest collision time or the smallest distance 
between scattering partners. In general, the potential impact that the screening of scattering partners by other test- 
particles has on physical scattering cross sections must be taken into account (see e.g. [64] and references therein). 
However, these effects can be avoided by requiring that the scattering partner selected for a given test-particle is always 
the closest of the potential candidates. While it is true that for matter in the hydrodynamic limit, the scattering 
cross section is very large and therefore the impact of screening might be negligible, we also find in our approach that 
the scattering partner determination according to the shortest distance proves itself as most suitable. We will discuss 
this in more detail in the next section. 



V. SIMULATION SETUP 



A main motivation of our work is to set up an algorithm for collision detection which is efficient and scales well 
to large simulations. The most computationally expensive step is the determination of scattering partners. A brute 
force approach would compare every particle to every other particle in the simulation. This is an 0(N 2 ) algorithm 
which could in theory run on P processors for an 0{N 2 / P) running time. An alternative is to assume that scattering 
partners are found based on nearest neighbors search. A simple approach to nearest neighbor is to use a sorting 
algorithm which is O(NlogN). Although there are methods for parallel sorting, this is still an expensive alternative 
and will not scale well to larger systems. For this work we choose to find scattering partners using spatial binning. 
Particles are placed into bins and then scattering partners are only searched for inside of a particular bin. The binning 
component of the algorithm takes linear time and will also scale linearly 0{N/P). The search for scattering partners 
is now determined by looking only at other particles in the same bin or close neighbors. If the maximum number of 
particles in a bin (iVbi n ) is much smaller than the total particles (iVbi n < logN < N) then this should significantly 
increase the computation time of the algorithm. Furthermore, this approach scales linearly for a total running time 
of 0(NN 2 JP). 

For this work the determination of scattering partners was parallelized using shared memory paralyzation (OpenMP). 
As is sketched in FigJI] (a) and (b), in a parallel setup, each CPU is assigned a bin and its neighboring bins. For a 
two-dimensional simulation the neighborhood consists of 8 cells, while in a 3D setup the number of neighboring bins 
is 27. Bins cannot be shared among processors since two CPUs might assign the same particle to different collision 
partners. The latter could lead to erroneous data in the particle's post-collision properties. Once a CPU has completed 
the collision partner search it is assigned a new bin. As is illustrated in Fig ljc) , there are two possibilities how to 
progress the parallel search. One option is to let all processors wait until the last CPU has finished its scattering 
detection. New bins are then assigned to all CPUs at the same time. The second possibility is that a processor can 
move independently to another bin once it has finished its scattering partner search. However, it must be ensured 
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FIG. 2: Flow chart of particle simulation 



that neither the new bin nor its neighborhood are part of the still ongoing collision detection of the other CPUs. For 
our work, we chose the first method. 

Figure [2] gives a rough outline of the code structure. The simulation starts with the initialization of the particle 
distribution. Hereby, particles are assigned their spatial coordinates, velocities, and further properties which are 
required in the simulation, such as mass. The simulation space has a given spatial extension. Due to the self-similar 
nature of shock wave problems neither the coordinates nor the velocities are given specific physical units. The binning 
grid is set up with a chosen number of bins in each dimension. We couple the time step size At to a particle velocity 
v and the width of a bin Aa; through: 

At 

At = 7?—. (12) 

v 

Hereby, v is a large particle velocity in the simulation and can either be coupled to the initial velocity distribution or 
be the maximal particle velocity at each time step. By setting the travelled distance to be only a small fraction t)< 1 
of the bin-width, we ensure that particles do not leave their scattering neighborhood within At. At the beginning of 
each time step, the particles are sorted into the bins, whereas, for each bin, we determine the corresponding number 
of particles iVbin- The latter is used to calculate the effective radius of each particle in two or three dimensions: 



Ax 2 Ax 3 f . 

rcff < 2D= 2AAW rcff ' 3D = V^AW ( ] 

Within a loop over all bins, the colli sion partner search is performed in parallel following the criteria for scattering 

Figure [3] shows the particle selection in one bin Abi n ~ O(10 4 ) for a chosen 



detection as discussed in section IV 



particle of interest according to the discussed collision tests: 
• Collision test 1: Check if <i m i n < (r e g a + ^eff.s) during At 

Determine the collision times i Cj 1 and i Ci 2 from a given A 



Collision test 2 



Collision test 3: Choose final collision partners according to their smallest distance 



With a chosen mean free path of A = O.OlAx, eq.(13) gives an effective particle radius of r G ff.2D = 5 x 10 _3 Aa; which 
explains why only a few particles in Fig. [3] pass the second collision test. Figure [5] gives a comparison between the 
total number of particles N in the simulation and the number of possible collisions. All simulations were performed 
for particles with the same initial velocity but random directions of motion. For a total particle number of N > 10 6 
we find that typically all particles in the simulation find a scattering partner. 

Once collision partners are assigned to each other, we perform the scattering in the center of mass frame of each 
colliding pair. This is currently done in a Newtonian formalism but will be extended to include effects of Special 
Relativity in the future, as have been done in previous works [SJ [7J [55] ■ The scattering angles can be calculated 
analytically for a given test-particle shape. However, due to the random nature of individual particle motion in bulk 
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FIG. 3: Three collision tests, as discussed in the text, performed in one bin for a particle number of Nbin ~ O(10 4 ), mean free 
path A = 0.01 As, velocity Vi = 300, and time step size dt = 0.1 Ax/vi. 



matter we chose a simpler approach by orienting the outgoing post-collision velocity vectors randomly and opposite 
to each other in space. 

For the scattering of final collision partners, different methods can be considered. Due to the large number of particles 
it is not unlikely that two particles, A and C, find the same collision partner B. Such situation can be resolved in 
different ways, for example: 

1. Scattering is only performed for two unambiguously determined collision partners 

2. Only one out of two conflicting scattering processes is performed 

3. Both scatterings are performed consecutively 

For the last case, the second scattering process involves the updated particle as collision partner, i.e. the latter 
has a new velocity vector as a result of the first interaction. Figure [4] shows collision partners with their pre- and 
post-collision velocity vectors for the above scenarios. The scattering partners are visualized through a connection by 
a thin dashed line. The mean free path was chosen to A = O.OlAx. With a number of particles per bin of N = 76, 
this leads to an effective radius of r c g ~ 0.66 Ax, which is large enough to ensure that typically every particle in the 
simulation should find a scattering partner. 

Note, that in the first scenario which is shown in figure [4] (a) many particles do not undergo collisions. In these cases 
an interaction partner was found. However, the search was in conflict with another particle's collision detection which 
determined the same particle as its scattering partner. As a result, no scattering was performed for all three particles. 
Such "ring structures" were noticed in previous studies and were shown to be able to lead to an underestimation of the 
collision rate |7j . Though Fig. [5] shows that the percentage of non-interacting particles seems to stays constant with 
growing particle number, the likeliness of conflicts in collision detection could increase for higher values of N. This 
could lead to problems in calculations aimed to describe dense systems in equilibrium. Instead of frequent particle 
interaction, the collision rate could decrease due to a growing number of scattering conflicts. 

Concerning the collision rate, the second approach for final collision partner search is more effective. Here, one of 
the two scatterings in conflict is performed and, as can be seen in Figj4] (b) and Figj5]the number of non-interacting 
particles is lower. Figure [4] (c) shows that for the third scenario we find that all particles are involved in scattering 
processes. With interactions in conflict being performed sequentially, the only possible scenario for which a possible 
interaction would be omitted is a situation in which particle A finds a scattering partner B, while the latter detects 
particle C as more suitable for a collision, and particle C does not interact at all. This case is very unlikely and we 
don't find its occurrence in our simulations. 

For all three scattering criteria, the final collision partners are chosen according to the smallest distance. For com- 
parison we also performed tests for the determination of the final scattering partner via the shortest collision time 
t c . The resulting collision partners are shown in Fig. [2] (d). It can be seen that, unlike for the smallest distance, 
interacting test-particles are generally farther away from each other. In case of large bin-widths or small values of A, 
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FIG. 4: Result of the scattering partner search with Nbin = 76, A = O.OlAa;, Vi = 300, and At — O.lAx/vi. Pre- and post- 
collision velocity vectors are illustrated by thick and thin lines, respectively. Collision partners are connected by a thin dashed 
line, (a) Unambiguously determined collision partners, (b) one out of two collisions in conflict is performed, (c) all collisions 
are performed, (d) Final collision partners are determined according to the shortest collision time. 
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FIG. 5: Possible collisions according to final collision criteria shown in Figs|4] (a) - (c) as a function of total particle number 
N. The number of final collisions is compared to all possible collisions in the simulation. The latter contain interactions which 
are in conflict with each other. 



such an algorithm for collision detection could decrease the spatial resolution in a simulation. For shock wave studies 
a failure to sufficiently spatially localize scattering partners could result in artificial broadening of the shock front. 
Therefore, to achieve a high resolution in our simulations we determine the final collision partner according to the 
closest distance. 

Once scattering has been performed for all colliding particles, we update the particle velocities and positions. For 
the latter, we distinguish between positive and negative collision times. As mentioned before, large particle densities 
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and small mean free paths can cause the particles' effective radii to overlap. This results in negative collision times 



according to eq.(10). In such a case we treat the particle collision as instantaneous and update the particle positions 
via: 

x(t + At) = x{t) + v ncw At, (14) 
where v ncw is the post-collision velocity. For particles with positive collision times t c the new position is: 

x(t + At) = x(t) + w id t c + v ncw (At - t c ) (15) 
where v id is the pre-collision velocity. 

In the next section we will present results of hydrodynamic shock wave studies performed with the kinetic test-particle 
code. 

VI. SIMULATIONS IN THE HYDRODYNAMIC LIMIT 

As mentioned in the introduction, an advantage of kinetic approaches is their ability to capture the behavior of 
matter for all Knudscn numbers. In this work we want to test the capability of our scattering partner search algorithm 
to reproduce hydrodynamic behavior, i.e. matter with small Knudsen numbers, especially with regard to the evolution 
of shock waves. The latter are often used as test cases to study the capability of a hydrodynamic code to capture 
steep gradients and discontinuities. Furthermore, most test problems possess analytic solutions which allow the user 
to quantify and compare the performance of different codes. The test suite which we present in this work consists of 
the following problems: the Sod test, a classical Riemann problem to show our code's ability to propagate shock waves 
and satisfy the shock jump conditions, and the Noh problem to test the conversion of kinetic in internal energy and 
shock wall-heating. We chose the Sod and Noh tests for the following reasons: the Sod test is generally considered 
to be an elementary benchmark test that all codes designed to model hydrodynamic Shockwave propagation should 
pass [55J and the Noh problem [57] tests the ability of our code to model the evolution of shocks in imploding systems 
similar to supernovae or ICF capsule systems that we aim to study. 

Our studies are typically carried out with ~ 10 7 test-particles placed in a simulation box which has reflective boundary 
conditions and is divided into (10 4 — 10 6 ) equally sized bins. To ensure hydrodynamic equilibrium by frequent 
interaction of particles, their mean free path is chosen to be a small fraction of the distance which the test-particles 
can traverse in a given timestep At. In the present study, interactions are only modeled in the form elastic scattering. 
As a consequence, the resulting matter behaves as a dense gas with the number of degrees of freedom / given by the 
number of dimensions in the simulation. The heat capacity ratio 7 = 1 + 1// becomes 7 = 2 for a two dimensional 
simulation and 7 ~ 1.67 in the case of three dimensions. The quantities which are calculated and compared to 
analytical solutions are the particle number density n, pressure p, and bulk velocity Vb- 

The density is given simply by the number of particles per corresponding volume Ny /V while the bulk velocity of 
matter is determined via: 

1 M 

V b = -j^-\J v l x + v t v + v l^ V bia = Y, V ^°" ( 16 ) 

^ i=l 

whereas a = x, y, and z. The particle velocities are also used in the determination of the pressure from the stress 
tensor of dense gases: [681169]: 

171 K<* - v b , a ) (Vi,p - v b .p) + ^ ^ mi n i- a Api <P ' ( 17 ) 



a/3 



with fij = fi — fj and Api — m (vi incw — Vi,o\d) being the change in momentum of particle i due to interaction with 



particle j. The stress tensor is composed of a kinetic and a potential part. Hereby, the first term in eq.(17l is the 
kinetic contribution. It dominates for gases and is induced by the momenta carried by the particles around their bulk 
motion. It is also the part of the stress tensor which results in the usual thermal pressure p = kbTN/V with kb being 
the Boltzmann constant. The second contribution is the potential part which arises due to momentum transfer from 



particle interaction and dominates the pressure in liquids. Applying eq.(17l to bulk matter with a volume V results 
in a corresponding pressure of: 

P2D = - V 2^ P3D = ~ V f • (18) 
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Velocity [v m ] Number of CPUs 

FIG. 6: Maxwell-Boltzmann distribution. M = 1.6 X 10 6 , 100 x 100 bins. From left top to right bottom: (a) A = 0.01 Ax (b) 
A = 0.1 Ax, (c) A = 1.0 Aa; (d) Scaling of the equilibration study with number of processors applying A = 0.01 Ax 

Given the symmetries of the simulated system, we calculate averages of n, p, or Vf, for a given interval of distance, 
e.g. distance Ax or radial distance Ar. In addition, number density, pressure, and bulk velocity are also determined 
as averages in sample bins. The grid used to define the sample bins does not have to be identical to the grid that 
defines the scattering bins. The number of sample bins can be chosen much larger for a high resolution data output 
if desired. 



A. Equilibration and CPU timing 

To study the equilibration rate in dependence of the test-particle mean free path, we set up a two-dimensional 
simulation space with N = 1.6 x 10 6 particles distributed homogeneously over 100 x 100 bins. All particles are 
initialized with the same initial absolute velocity v- ln and random orientation of the velocity vector. Varying the mean 
free path from A = O.OlAir, O.lAa;, to l.OAai, we determine the particle velocity distribution over the first ten time 
steps. For systems in equilibrium the latter should promptly reach the Maxwell-Boltzmann distribution. Figures 
[6[a)-[6jc) show the number of particles as a function of absolute velocity for different values of A. It can be seen that 
simulations with A = 0.01 Ax and A = 0.1 Ax obtain the Maxwell Boltzmann velocity distribution within the first 
6-10 time steps, while for mean free paths A > 1.0 Ax the equilibration time is longer. This is a motivation for us to 
apply small values of A < 0.1 Ax in simulations which aim to describe equilibrium configurations. 
As discussed in sections |IV| and [Vj our scattering partner search algorithm is designed to work in parallel. To test 
the speed-up we performed timing tests for the first ten timesteps of the above velocity equilibration study using 
A = 0.01 Aa;. The computational speed-up is defined as S = T#cpu /Ti, where T#cpu is the wall-clock performance 
time of an algorithm with #CPU number of processors, while T\ is the time when only one processor is used. The 
results for #CPU = 1 — 32 are shown in figure [6jd) which gives as comparison also the case of ideal speed-up 
with S = #CPU. Though the performed simulations show a deviation from the latter with increasing number of 
processors the overall speed-up is still close to an ideal behavior and therefore promising for implementations of our 
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x bins x bins x bins 

FIG. 7: 2D Sod shock test: TV = 2.4x 10 7 test-particles are distributed over 400x400 bins with A = O.OlAa; and At = 0.25Ax/vr. 
(a) Particle number per bin, (b) pressure, and (c) bulk velocity with developed shock profiles at timestep t — 350 At 

algorithm on larger systems with thousands of cores. It should however be noted that the CPU scaling shows such 
a good performance due to the homogeneous distribution of test-particles and thereby equal load of all processors. 
For physical problems which involve a large number of particles concentrated in only a few bins and therefore not 
distributed over all CPUs the scaling will be worse. Possible solutions to ensure a more homogeneous processor load 
include the usage of adaptive grids or the individual CPU update of the active bins as discussed in section [V] 



B. Sod test 



We start our series of shock wave studies with the Sod shock-tube test [66]. It is a classical a Riemann problem 
and has been performed by various hydrodynamic and kinetic algorithms (see e.g. |70H77| ) . The simulation space is 
divided into two equally sized parts with different pressures and densities on both sides. The ratios of the left- and 
right-hand side pressures, number densities, and bulk velocities are set to: 

{n,p,v b ) L = (1,1,0), (n,p,v b ) R = (0.125,0.1,0). (19) 

Initially, both box partitions are separated. At the start of the simulation, the separation wall is removed which leads 
to the development of a shock front that propagates from the high density into the low density matter. The shock 
front is followed by a contact discontinuity and a rarefaction fan. The latter is visible in the pressure as well as the 
density profiles, while the contact discontinuity manifests itself only in the density distribution. 

Our first simulation is carried out with N = 2.4 x 10 7 particles in two dimensions. The size of the simulation space 
spans over < x,y < 7 and is divided into 400 x 400 bins. We initialize the particle velocities with = 300 and 



vl = 268.32, while the particle density ratios are set up as in eq.(19). With a chosen particle mean free path of 
A = 0.01 Ax and a timestep size of At = 0.25 Ax/vr the particle velocities equilibrate to the Maxwell-Boltzmann 
distribution within the first time steps. After 20 time steps the partition wall is removed and we observe the evolution 
of a shock front which propagates from the right into the low density and low pressure region on the left. 
Figure [7] shows the particle number per bin, the pressure, and the bulk velocity at timestep t = 350 At. To compute 
the corresponding profiles as a function of distance x, we average the density and pressure over bins in the y-direction, 
while the bulk velocity is computed via cq.(16) with averaged v x and v y . The resulting profiles are then compared 



to the analytic solutions which are obtained using the exact-riemann.f code |78| . We plot the results for density, 
pressure, and bulk velocity in Figs.[8ja) - [8jc) , whereas the density and pressure are normalized by their initial values 
riL and p^. We find that the simulation agrees well with the analytic solutions. Disagreement can be found in the 
sharpness of the contact discontinuity in the particle number density. The latter is accompanied by an increase in the 
bulk velocity, an effect which has also been observed in other kinetic approaches [72l [79] and could be connected to 
the particle mean free path 80J. Statistical fluctuations can be found in all quantities but are expected to decrease 
for larger particle numbers. 



Results for the three-dimensional Sod shock tube test are shown in Fig. [9] and Fig. 10 For better comparison with the 



two dimensional calculation, we set up the same particle densities by dividing the simulation space into 400 x 20 x 20 
bins and keeping the bin-widths as well as the particle numbers the same as in the 2D case. The resulting density, 



pressure, and bulk velocity profiles are presented in Figs. |10[a) - 10 c). The results are qualitatively very similar 



to the two dimensional simulations whereas the heights of the pressure and density plateaus are modified due to 
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FIG. 8: 2D Sod shock test as in Fig. [7] Profiles of (a) normalized density, (b) normalized pressure, and (c) bulk velocity at 
t — 350 At are shown together with the analytical solutions. 



(a) Particles per bin (b) Pressure [10 9 ] (c) Bulk velocity 




FIG. 9: 3D Sod shock test: N = 2.4 x 10 7 test-particles are distributed over 400 x 20 x 20 bins with A = 0.01 Ax and 
At — 0.25 Ax/vr. (a) Particle number per bin, (b) pressure, and (c) bulk velocity with developed shock profiles at timestep 
t = 350 At. 




FIG. 10: 3D Sod shock test as in Fig. [9] Profiles of (a) normalized density, (b) normalized pressure, and (c) bulk velocity at 
t — 350 At are shown together with the analytical solutions. 
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FIG. 11: Planar Noli test in 2D: N = 3.2 x 10 7 test-particles distributed over 300 x 150 bins with A = 0.01 Aa; and At = 
0.25Ax/ui n . (a) Particle number per bin, (b) pressure, and (c) bulk velocity at timestep t — 380 At with developed shock 
profiles. 



the different heat capacity ratio of 7 ~ 1.67 for a three dimensional system. Again, it can be seen in Fig. flO] that 
the analytical solutions are well reproduced by the test-particle approach. This has also been demonstrated in other 
kinetic simulations |72H77j . With that, we turn to the Noh test which is a more challenging problem for hydrodynamic 
simulations. 



C. Noh test 



The Noh or Newtonian wall shock problem |67j is an important test for approaches which aim to study collapsing, 
i.e. imploding systems. While being performed by hydrodynamic codes [81] , the Noh problem has not been analyzed 
in many particle based approaches, and is therefore of special interest for our study. The test is initialized in form 
of a gas with homogeneous density nj n that streams with a uniform velocity vi n towards a wall of the simulation box 
or its origin. As matter starts to pile up, it converts kinetic into internal energy. An accretion shock front forms 
enclosing the shocked matter, whereas its radius as a function of time is given by: 

r s hock(*) = ^(7-l)«mi- (20) 
The resulting density profile can be written as a function of distance r: 

/7+lV 

n(r) = n , r < r shock (21) 



7-1 

n(r) = n ( 1 + — J , r > r shock . (22) 

Hereby, d gives the geometry of the system with d — 1 for a planar, d — 2 for a cylindrical, and d = 3 for a spherical 
setup. Many hydrodynamic simulation codes experience problems with the Noh test due to anomalous wall-heating 
which decreases the density of matter close to the wall of the simulation box [57J [5T] ■ 

We performed the Noh test in planar as well as cylindrical geometry with a total particle number of N = 3.2 x 10 7 
and Wi n = 300. For planar simulation we divide the simulation space into 300 x 150 bins with < x < 5.26 and 
< y < 2.625. The mean free path is chosen to A = 0.01 Ax with a timestep size At = 0.25 Ax/fj n . The number 



of particles, pressure, and bulk velocity per bin at t = 300 At are given in Fig. 11 The figures show the shocked 



high density matter composed of particles which are streaming from the right-hand side of the simulation box and 



accumulate along the left wall. Figures 12 'a) -|l2[c) give the corresponding normalized particle densities, pressures, 
and bulk velocities averaged over the y-dimension. Generally the results of the simulation are in good agreement with 
the analytical solution. The latter was obtained by the noh.f code from [75]. The particle density of the shocked 
matter is slightly under-predicted and the shock front is again broader than in the analytic solution. However, as can 
be seen in Fig. |13[ the density profiles are sensitive to Abj n and approach the analytic solution when its values are 
increased. 
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FIG. 12: Planar Noli test in 2D as in Fig. 11 Profiles of (a) normalized density, (b) normalized pressure, and (c) bulk velocity 



at t = 380 At are shown together with the analytical solutions. 
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FIG. 13: Noh test for a two dimensional setup as described in Fig. [T2] with different values of particles per bin Abin 
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FIG. 14: Cylindrical Noh test in 2D: N = 3.2 x 10 7 test-particles distributed over 200 x 200 bins with A = 0.01 Ax and 
At — 0.25Ai/wi n . (a) Particle number per bin, (b) pressure, and (c) bulk velocity at timestep t = 270 At with developed shock 
profiles. 



15 



(a) Particle density n/ri;, 



10 

9 
8 
7 
6 
5 
4 
3 
2 
1 



"v*w^«~-=r| Analytic 



Simulation 



0.2 0.4 0.6 0.8 1 1.2 1.4 
r 



1.2 

0.8 
0.6 
0.4 
0.2 


-0.2 



(b) Pressure p/p ;i 



(c) Bulk velocity Vj/Vy 



Simulation 
Analytic 



0.2 0.4 0.6 0.8 1 1.2 1.4 
r 



1.2 
1 

0.8 
0.6 
0.4 
0.2 






Simulation 




Analytic 


- 




, * 





0.2 0.4 0.6 0.8 1 1.2 1.4 
r 



FIG. 15: Cylindrical Noh test in 2D as in Fig. 14 Profiles of (a) normalized density, (b) normalized pressure, and (c) bulk 



velocity at t = 270 At are shown together with the analytical solutions. 



In Fig. 14 and 15 the results of the cylindrical Noh test are shown. The simulations space is now divided into 200 x 200 
bins with < x,y < 3.5, whereas timestep size and mean free path are the same as for d = 1. The particle number 
per bin, pressure, and bulk velocity are shown in Fig. |14| Different to the Noh test in planar geometry, particles 
now stream with uniform radial velocity towards the origin of the simulation. The density, pressure, and velocity 
averages are also taken over the radial distance r. Again, we find good agreement between analytical prediction and 
kinetic simulation as shown in Fig. [TB^a) - Fig. [l5jc), while small deviations are present in the particle density and 
the width of the shock front. In addition, for the shocked matter, the bulk velocity in Fig. [l5|c) seems to be above 
zero. The reason for the this lies in our approach to calculate Vj,- For the cylindrical Noh test we first determine 
the bulk velocity according to eq.(16l for every bin and then calculate its radial average. As a consequence, possible 
negative and positive fluctuations in v x and v y for different bins do not average out but are summed and thereby lead 
to positive fluctuations of Vb- The latter, averaged over radial distance, results in the observed bulk velocity above 
zero. An increase of the bin-width or the number of test-particles per bin A^m should improve the statistics and 
reduce the bulk velocity. An alternative approach is to determine the radial velocity: 



1 



M 



i—l 



(23) 



instead. It is a function of radial distance r and can have positive and negative values. Thereby fluctuations of v r 
can cancel out when averaged over several bins, leading to the expected v r ~ for shocked matter. However, for 
both, the planar as well as cylindrical Noh tests, it is noteworthy that no signs of wall-heating are observed. It will be 
interesting to see if this observation stays invariable for even higher particle numbers and how it compares to other 
particle based simulations. 



VII. SIMULATIONS OF NON-EQUILIBRIUM SYSTEMS 



As mentioned in the introduction, an important advantage of kinetic schemes is their ability to simulate systems 
which have large Knudsen numbers and are therefore not equilibrium. In this section, we want to test this regime 
within our approach using particle mean free paths that are larger than a characteristic length scale of the system 
A > L. We choose L = 3 Ax since it defines a particle's scattering neighborhood and study two representative values 
of A = 4 Ax and A = 32 Ax. Both are applied to shock tests which were previously discussed in sections |VIB| and 
|VI Cj whereas we compare the results to the equilibrium regime with A = 0.01 Aa; and the free-streaming case. For 
the latter the particles do not interact with each other and their mean free path can be assumed as infinitely large. 
As representative cases for the shock simulations we choose the 2D Sod test and the 2D planar Noh test, whereas we 
expect the results for the 3D Sod test and the cylindrical Noh test to be qualitatively similar. Mean free path studies 
of the Sod test have also been performed by e.g. [SJ [50] whereas we are not aware of any similar studies with the Noh 
problem. 
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(a) Particle density n/n L (b) Pressure p/p L (c) Bulk velocity v b 




FIG. 16: 2D Sod shock test as in Fig. [7] Profiles of (a) density, (b) pressure, and (c) bulk velocity at t — 350 At. The mean 
free path is varied from A = 0.01 Aa;, 4 Aa; to 32 Ax. The simulation using free streaming particles is marked by Fs. 
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FIG. 17: 2D Noh shock test as in Fig. 12 in planar geometry and N = 1.6 x 10 7 . Profiles of (a) density, (b) pressure, and 
(c) bulk velocity at t — 200 At. The mean free path is varied from A = 0.01 Aa;, 4 Aa; to 32 Aa;. The simulation using free 
streaming particles is marked by Fs. 



We perform the Sod test with N = 2.4 x 10 7 test-particles, 400 x 400 bins, and —3.5 < x,y < 3.5. The initialization 
process is the same as in section |VI B| For the first twenty timesteps the particle system equilibrates to the corre- 
sponding Maxwell-Boltzmann velocity distributions whereas the test-particle mean free paths are A = 0.01 Ax. At 
t = 20 At, the partition wall between the high and low density matter is released and the particle mean free paths 
are set to their new values of A = 4 Aa; and A = 32 Ax, respectively. The results are shown in Fig. 16 'a) - I 
We plot the density, pressure, and bulk velocity averaged over the y-dimension. It can be seen that for A = 4 A 



16jc). 
x the 



density, pressure, and velocity profiles are significantly washed out, however, they still indicate the general features 
of the Sod shock. For A = 32 Ax, on the other hand, almost no shock features are visible. The density profile is very 
similar to the free-streaming regime, whereas the pressure profiles still reproduces a small plateau near the contact 
discontinuity. 

Similar behavior can be observed for the Noh test in Fig. |l7{a) - Fig. 17 c). As before, we plot the density, pres- 
sure, and bulk velocity profiles as functions of the distance x. The simulation space is set up with N = 1.6 x 10 7 
test-particles in 300 x 150 bins, where < x < 5.25 and < y < 2.625. The timestep size and the initial particle 



velocity are the same as in section VIC Like in the Sod mean free path test, it can be seen that with increasing 
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values of A, the shock front broadens. This can be attributed to particles which, once they reach the shocked matter 
region, are not trapped by in-flowing particles, but can escape due to their large mean free path. While the peak 
density and pressure for A = 4 Ax are higher than in the equilibrium case, the opposite behavior can be found for 
A = 32 Ax. A possible reason is that for A = 4 Ax more particles reach the wall of the simulation box at x = 0. 
Here, they can become trapped due to scattering with in-flowing particles which leads to a higher particle number. 
If the mean free path is set to larger values of A = 32 Ax, the majority of particles does not interact and, after being 
reflected on the walls of the simulation box, propagates back to larger distances. As can be expected, this behavior 
is even more pronounced in the free streaming regime. This results in densities twice as large as the initial value and 
zero bulk velocity. The latter is a consequence of the cancellation between the bulk velocities of in- and out-flowing 
test-particles. The non-zero pressure stems from the first term in eq.(17) and is due to the non- vanishing individual 
particle motion versus the zero total bulk velocity. 



VIII. SUMMARY AND CONCLUSION 



In this work we introduce and discuss a Monte Carlo kinetic scheme which simulates the phase space evolution of 
systems that can move in and out of hydrodynamic equilibrium. The code operates with test-particles that interact 
with each other via two-body collisions with variable particle mean free paths. To achieve high computational 
efficiency, the code is written in a parallel form and applies spatial binning of test-particlcs. We determine the 
interaction partners using the distance of closest approach while the final collision partners are chosen according to 
the shortest distance between them. Both deviations from the usual Direct Simulation Monte Carlo techniques, which 
determine interaction partners randomly from a scattering cell, are aimed to enhance spatial accuracy which could 
otherwise be compromised due to scattering partners which are far away from each other. 

For first performance studies of our kinetic scheme we apply timing and mean free path tests with a simple gas of test- 
particles moving in two dimensions and interacting with each other via binary collisions. With equal initial absolute 
values, the test-particle velocities quickly equilibrate to the Maxwell-Boltzmann distribution when the particle mean 
free path is chosen small with respect to the distance which can be travelled by a particle within one timestep. Larger 
mean free paths on the other hand, lead to longer equilibration times. In addition, for this simple setup we find 
that the computational speed-up of the scattering partner search with the number of processors shows almost ideal 
behavior. 

To further explore the capability of our kinetic scheme to reproduce hydrodynamic behavior as well as handle steep 
gradients and discontinuities, we apply shock wave tests in fixed-volume and imploding systems. These are the Sod 
and the Noh tests which we perform with 2.4 x 10 7 and 3.2 x 10 7 test particles, respectively, and compare the evolution 
of the resulting shock structures to the corresponding analytic solutions. Generally, we find good agreement between 
the predicted shock profiles and the results of the kinetic simulations. Common deviations are found in the sharpness 
of the shock fronts whereas these improve when larger particle numbers are used. Finally, we test the ability of the 
kinetic particle scheme to evolve systems which are not in hydrodynamic equilibrium. We perform the Sod and the 
Noh tests applying particle mean free paths which are much larger than the distance which a test-particle can travel 
within one timestep. As expected, we find that larger particle mean free paths lead to a more pronounced broadening 
of the shock front whereas the solutions approach the free streaming limit. 

We can conclude that the particle kinetic scheme is capable to reproduce hydrodynamic behavior, especially with 
regard to shock wave evolution. More importantly, it can also simulate systems which are not in hydrodynamic 
equilibrium by applying large particle mean free paths. With this, it is an interesting tool to study systems such 
as core-collapse supernovae where neutrinos traverse from the trapped to the free streaming regime as well as the 
dynamics of ICF capsules where components of the nuclear fuel, such as thermal deuterons, posses mean free paths 
which are far too large relative to characteristic length scales of the region of thermonuclear ignition to assume 
hydrodynamic equilibrium. 



IX. OUTLOOK 



Future applications of our code in e.g. core-collapse supernova simulations will require further performance studies 
such as additional shock tests in two and three dimensions as well as studies concerning the ability of the particle 
approach to reproduce turbulent fluid motion, for example Rayleigh- Taylor and the Helmholtz instabilities. The search 
for scattering partners has to be upgraded to be able to handle particles moving at relativistic speeds. Also, three 
dimensional simulations will require higher particle numbers than were applied in the current study of N ~ 10 8 — 10 9 
test-particles. To handle such large particle numbers, the scattering partner search has to be set up with distributed 
memory parallelization to enable the usage of S> 10 2 processors. Corresponding CPU timing tests have to be performed 
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whereas large simulations would benefit from additional extensions such as adaptive grid and time step sizes for a 
better handling of shocks and large particle numbers. Another important implementation for future applications 
are the usage of additional test-particle interaction forces and realistic particle scattering cross-sections, e.g. nuclear 
forces and neutrino interactions with nuclear matter in core-collapse supernovae, as well as the handling of long-range 
forces such as gravitation. Finally, we plan to provide the scattering partner search code in form of a library to be 
applicable in a variety of studies of different hydrodynamic and transport problems. 
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